Hands-on Ex 3

Network Constrained Spatial Point Pattern Analysis

Author

Stephen Tay

Published

September 8, 2024

Modified

September 8, 2024

1. Overview

Network-constrained spatial point pattern analysis (NetSPAA) focuses on studying events that occur along networks, such as roads and rivers, rather than across open spaces. Unlike traditional spatial analysis, which assumes events can happen anywhere and measures distances in straight lines, this method addresses phenomena restricted to networks, like traffic accidents or store locations. It is especially important for its accurate in detecting event clusters that occur along these networks.

The spNetwork package is a crucial geospatial R tool for conducting NetSPAA, specifically for:

  • Network Kernel Density Estimation (NKDE), and
  • Network-based G-function and K-function analysis.
pacman::p_load(sf, spNetwork, tmap, tidyverse)

2. Importing & Transforming Data

This study analyses the spatial distribution of childcare centers in the Punggol Planning Area. We will import and work with two geospatial datasets:

  • Punggol_St: A line feature dataset representing the road network within the Punggol Planning Area.
  • Punggol_CC: A point feature dataset representing the locations of childcare centers within the Punggol Planning Area.

After importing the road network data, it is essential to review the dataframe, verify the coordinate reference system (CRS) is correct, and check for any duplicated geometries.

network <- st_read(dsn="data/geospatial", 
                   layer="Punggol_St")
Reading layer `Punggol_St' from data source 
  `/Users/stephentay/stephentay/ISSS626-Geospatial-Analytics/Hands-on_Ex/Hands-on_Ex03/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TM
st_crs(network)
Coordinate Reference System:
  User input: SVY21 / Singapore TM 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
network
Simple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TM
First 10 features:
     LINK_ID                   ST_NAME                       geometry
1  116130894                PUNGGOL RD LINESTRING (36546.89 44574....
2  116130897 PONGGOL TWENTY-FOURTH AVE LINESTRING (36546.89 44574....
3  116130901   PONGGOL SEVENTEENTH AVE LINESTRING (36012.73 44154....
4  116130902   PONGGOL SEVENTEENTH AVE LINESTRING (36062.81 44197....
5  116130907           PUNGGOL CENTRAL LINESTRING (36131.85 42755....
6  116130908                PUNGGOL RD LINESTRING (36112.93 42752....
7  116130909           PUNGGOL CENTRAL LINESTRING (36127.4 42744.5...
8  116130910               PUNGGOL FLD LINESTRING (35994.98 42428....
9  116130911               PUNGGOL FLD LINESTRING (35984.97 42407....
10 116130912            EDGEFIELD PLNS LINESTRING (36200.87 42219....
any(duplicated(network))
[1] FALSE

After importing the point data, it is essential to review the dataframe, verify the CRS is correct, and check for any duplicated geometries. Since the point geometries are in 3-dimensional coordinates, they must be transformed to 2-dimensional coordinates using the st_zm() method.

childcare <- st_read(dsn="data/geospatial",
                     layer="Punggol_CC")
Reading layer `Punggol_CC' from data source 
  `/Users/stephentay/stephentay/ISSS626-Geospatial-Analytics/Hands-on_Ex/Hands-on_Ex03/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range:       zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM
st_crs(childcare)
Coordinate Reference System:
  User input: SVY21 / Singapore TM 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
childcare
Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range:       zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM
First 10 features:
      Name                      geometry
1   kml_10 POINT Z (36173.81 42550.33 0)
2   kml_99 POINT Z (36479.56 42405.21 0)
3  kml_100 POINT Z (36618.72 41989.13 0)
4  kml_101 POINT Z (36285.37 42261.42 0)
5  kml_122  POINT Z (35414.54 42625.1 0)
6  kml_161 POINT Z (36545.16 42580.09 0)
7  kml_172 POINT Z (35289.44 44083.57 0)
8  kml_188 POINT Z (36520.56 42844.74 0)
9  kml_205  POINT Z (36924.01 41503.6 0)
10 kml_222 POINT Z (37141.76 42326.36 0)
any(duplicated(childcare))
[1] FALSE
childcare$geometry <- st_zm(childcare$geometry)
childcare
Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
Projected CRS: SVY21 / Singapore TM
First 10 features:
      Name                  geometry
1   kml_10 POINT (36173.81 42550.33)
2   kml_99 POINT (36479.56 42405.21)
3  kml_100 POINT (36618.72 41989.13)
4  kml_101 POINT (36285.37 42261.42)
5  kml_122  POINT (35414.54 42625.1)
6  kml_161 POINT (36545.16 42580.09)
7  kml_172 POINT (35289.44 44083.57)
8  kml_188 POINT (36520.56 42844.74)
9  kml_205  POINT (36924.01 41503.6)
10 kml_222 POINT (37141.76 42326.36)

3. GeoVisualisation of Childcare Centers

Before starting the analysis, it is good practice to visualize the geospatial data.

plot(st_geometry(network))
plot(childcare,add=T,col='red',pch = 19)

tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(childcare) + 
  tm_dots() + 
  tm_shape(network) +
  tm_lines()
tmap_mode('plot')
tmap mode set to plotting

4. Network Kernel Density Estimation (NKDE) Analysis

4.1 Prepare lixels

Before computing NKDE, the SpatialLines object must be segmented into lixels with a specified minimum distance. The lixel length (lx_length) is set to 700m, and the minimum length (mindist) is set to 350m. After segmentation, if the final lixel is shorter than the minimum distance, it is merged with the previous lixel.

lixels <- lixelize_lines(network, 
                         lx_length = 700, 
                         mindist = 350)

4.2 Compute lines’ center point

Next, lines_center() from the spNetwork package will be used to generate a SpatialPointsDataFrame with points positioned at the center of each line, based on the line’s length, as shown in the code below.

samples <- lines_center(lixels) 

4.3 Perform NKDE

We are ready to computer the NKDE. There are 3 methods in computing NKDE:

The simple method, proposed by Xie and Yan (2008), differs from classical KDE in two key ways:

  • Events are snapped to a network.
  • Distances between sampling points and events are measured along the network (i.e., not using Euclidean distances).

Let’s examine the NKDE plot using the simple method below!

simple_densities <- nkde(lines = network, events = childcare,
                  w = rep(1, nrow(childcare)), samples = samples,
                  kernel_name = "quartic", bw = 300, 
                  div= "bw", method = "simple", 
                  digits = 1, tol = 1,
                  grid_shape = c(1,1), max_depth = 8,
                  agg = 5, sparse = TRUE,
                  verbose = FALSE)
samples$simple_density <- simple_densities*1000
lixels$simple_density <- simple_densities*1000
tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(lixels) +
  tm_lines(col="simple_density") +
  tm_shape(childcare) +
  tm_dots()
tmap_mode('plot')
tmap mode set to plotting

The simple method has a limitation: the event’s mass is multiplied by the number of edges at an intersection, causing the kernel not to integrate to 1. To address this, Okabe and Sugihara (2012) proposed the discontinuous NKDE, which was further extended by Sugihara, Satoh, and Okabe (2010) for networks with cycles shorter than the bandwidth.

See the NKDE plot using the discontinuous method below!

discontinuous_densities <- nkde(lines = network, events = childcare,
                  w = rep(1, nrow(childcare)), samples = samples,
                  kernel_name = "quartic", bw = 300, 
                  div= "bw", method = "discontinuous", 
                  digits = 1, tol = 1,
                  grid_shape = c(1,1), max_depth = 8,
                  agg = 5, sparse = TRUE,
                  verbose = FALSE)
samples$discontinuous_density <- discontinuous_densities*1000
lixels$discontinuous_density <- discontinuous_densities*1000
tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(lixels) +
  tm_lines(col="discontinuous_density") +
  tm_shape(childcare) +
  tm_dots()
tmap_mode('plot')
tmap mode set to plotting

Discontinuous NKDE can be counterintuitive, resulting in sharp differences between density values across the network, particularly problematic in networks with many closely spaced intersections. The continuous method adjusts NKDE values at intersections to ensure the kernel integrates to one, applying a backward correction to maintain continuity in density values. However, this method is recursive and computationally time-consuming.

See the NKDE plot using the continuous method below!

continuous_densities <- nkde(lines = network, events = childcare,
                  w = rep(1, nrow(childcare)), samples = samples,
                  kernel_name = "quartic", bw = 300, 
                  div= "bw", method = "continuous", 
                  digits = 1, tol = 1,
                  grid_shape = c(1,1), max_depth = 8,
                  agg = 5, sparse = TRUE,
                  verbose = FALSE)
samples$continuous_density <- continuous_densities*1000
lixels$continuous_density <- continuous_densities*1000
tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(lixels) +
  tm_lines(col="continuous_density") +
  tm_shape(childcare) +
  tm_dots()
tmap_mode('plot')
tmap mode set to plotting

5. Network-constrained G- and K-function analysis

A test for complete spatial randomness (CSR) can be performed using the network-constrained G- and K-functions.

The test hypotheses are as follows:

  • H0: The distribution of childcare services in Punggol are randomly distributed along the road network.
  • H1: The distribution of childcare services in Punggol are not randomly distributed along the road network.

Alpha is set at 0.05, with 1000 Monte-carlo simulations.

set.seed(2024)
kfun_childcare <- kfunctions(network, 
                             childcare,
                             start = 0, 
                             end = 1000, 
                             step = 50, 
                             width = 50, 
                             nsim = 1000, 
                             resolution = 50,
                             verbose = FALSE, 
                             conf_int = 0.05)

The blue line represents the empirical network K-function of the childcare centers in Punggol, while the gray area shows the results of 1,000 simulations within the 2.5% - 97.5% confidence interval. Since the blue line falls below the lower boundary of the envelope between the 200-400m distance range, it indicates a regular pattern of childcare center locations along the Punggol road network.

kfun_childcare$plotk

The blue line represents the empirical network G-function of the childcare centers in Punggol. Since the blue line falls below the lower boundary of the envelope in the 100-200m distance range, it suggests a regular pattern of childcare center locations along the Punggol road network.

kfun_childcare$plotg

6. Summary

Network-constrained spatial point pattern analysis revealed a regular distribution of childcare centers along the Punggol road network.